Interpolating compact binary waveforms using the singular value decomposition 
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Compact binary systems with total masses between tens and hundreds of solar masses will produce 
gravitational waves during their merger phase that are detectable by second-generation ground- 
based gravitational- wave detectors. In order to model the gravitational waveform of the merger 
epoch of compact binary coalescence, the full Einstein equations must be solved numerically for the 
entire mass and spin parameter space. However, this is computationally expensive. Several models 
have been proposed to interpolate the results of numerical relativity simulations. In this paper 
we propose a numerical interpolation scheme that stems from the singular value decomposition. 
This algorithm shows promise in allowing one to construct arbitrary waveforms within a certain 
parameter space given a sufficient density of numerical simulations covering the same parameter 
space. We also investigate how similar approaches could be used to interpolate waveforms in the 
context of parameter estimation. 



I. INTRODUCTION 

Searches for gravitational waves from binary black 
holes with total masses between tens and hundreds of 
solar masses benefit from the complete model of the 
gravitational waveform obtained by numerical relativ- 
ity [ , ]. Numerically solving Einstein's equations is 
now quite routine [3-9], yet still computationally bur- 
densome. Reference [ ] suggests that there is a finite 
density of numerical simulations that would adequately 
cover the parameter space for certain ground-based de- 
tectors. In this work we explore this concept and extend 
the numerical techniques presented in [ ] and [ ], to 
interpolation of template waveforms using the singular 
value decomposition. This should allow for the construc- 
tion of gravitational waveforms with parameters between 
the numerically generated waveforms. 

The idea of interpolating gravitational waveforms has 
existed for over a decade. Interpolation of waveforms 
generated by post-Newtonian techniques was described 
in [13] and [14]. In these references analytic formulae for 
waveform interpolation were derived for particular PN 
models. Since 2005 the numerical relativity community 
has been generating a substantial number of gravitational 
waveforms for the coalescence of binary black holes [3-9] . 
Interpolation of these waveforms has been accomplished 
primarily by (i) phenomenologically fitting the simula- 
tions to closed- form expressions [15-17] or (ii) by numer- 
ically solving simpler differential equations that capture 
the orbital dynamics combined with numerical stitching 
of the ringdown phase [18-23]. In this work we propose a 
different approach to interpolate a set of template wave- 
forms. This approach does not require careful tuning 
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of fitting formulae or stitching of waveforms and can be 
applied to any waveform set of sufficient density. 

This paper is organized as follows. First, we describe 
the technique for interpolating waveforms via the singu- 
lar value decomposition. Second we apply the technique 
to a set of waveforms containing all phases of the compact 
binary coalescence, inspiral, merger and ringdown. Third 
we discuss how these results might be applied to the con- 
duction of waveform families, ongoing gravitational wave 
searches, and parameter estimation. 



II. INTERPOLATION TECHNIQUE 

It was shown in [ ] that the singular value 
decomposition (SVD) reduces the number of template 
waveforms needed to search a given parameter space. Ad- 
ditionally, [ ] showed that arbitrary waveforms within 
the parameter space could be reconstructed from the 
SVD of a sufficiently dense template bank. Here we 
demonstrate a method to directly obtain approximate 
reconstruction coefficients for arbitrary waveforms in the 
parameter space via interpolation. Consider a wave- 
form family h(x, y) described by the physical parameters 
(x, y) : and consider a set of these waveforms enumerated 
by the index a drawn from a region of the parameter 
space, h(x a ,y a ). Recall that a SVD of these waveforms 
allows each to be written as a linear combination of basis 
waveforms u M with coefficients M^(x ai y a ) 

H x <*,y*) = ^M^x^y^u^, (1) 

where, in the formalism of [ ] and [ 2], M^{x ai y a ) := 
cr /i(^(2o;-i)/i + iv {2ol)ll) is the ath combination of singu- 
lar values cr^ and reconstruction coefficients v^a-i)^ an d 
v (2a)fi- Recall also that waveforms with arbitrary physi- 
cal parameters from the same region of parameter space 
can also be reconstructed using the basis vectors u M by 
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projecting the waveforms onto the basis vectors to ob- 
tain the reconstruction coefficients — a computationally 
expensive procedure, 



h(x, y) « ^2(h(x, y) • u M ) u M . 



(2) 



This can be used to define the arbitrary reconstruction 
coefficients as 



M^{x,y) h(x,y) 



(3) 



We seek the set of interpolated reconstruction coefficients 
M'^x^y) that can approximately reconstruct an arbi- 
trary waveform from that region of parameter space. 

Compact binary gravitational waveforms with negligi- 
ble effects from spin and eccentricity are characterized by 
their component masses. We will assume for concreteness 
a two parameter family of waveforms h(x, y) where x and 
y are M and q, respectively, where M = mi + rri2 is the 
total mass of the system and q = rai/7712 is the mass 
ratio of the system. 

Chebyshev polynomials of the first kind are known to 
be good for interpolation, however other interpolation 
schemes are also possible. We start with a set of basis 
vectors Uj covering the desired region of parameter space. 
We choose a net of points, scaled such that each dimen- 
sion covers the interval [—1,1], located at the J ma xth 
order Chebyshev nodes. For a single dimension, these 
nodes occur at the locations 



Xj = COS I 7T- 

Jrr 



1 



(4) 



where j ranges from to J ma x- This choice of net re- 
duces Runge's phenomenon when used with the Cheby- 
shev polynomials, which, for a single dimension, are given 
as 



Tj(x) 



(x - Vx 2 - 1) J + (x + Vx 2 - 1) J 



2w 



(5) 



where w = + Sjo)(J max + l)/2 is a normalization 
factor for the polynomials and Sjo is the Kroenecker 
delta. Both xj and w depend on the choice of J ma x 7 
however for ease of notation we will leave this implied. 
The polynomials Tj(x) satisfy the discrete orthogonality 
condition 



J max 

£ 

3=0 



Tjix^Tjixj) = 6. 



(6) 



It is straightforward to extend this to higher dimensions. 

In order to obtain the reconstruction coefficients for 
these locations, we project waveforms from these loca- 
tions onto the basis vectors. From the values on this 
net, we interpolate to other positions in parameter space 
using 2D-Chebyshev interpolation for each set of recon- 
struction coefficients M fl (x^y). Specifically, these values 



are projected onto the Chebyshev polynomials 

^max -^max 

(x k ,yi). (7) 



k=0 1=0 



This results in coefficients for the 2D-Chebyshev poly- 
nomials which can be used to evaluate the interpolated 
reconstruction coefficients M^(x, y) at other points in pa- 
rameter space 



K&y) =EE C K L»T K (x)T L (y). (8) 



k=o L=0 



In the next section we explore this approximation tech- 
nique using gravitational waveforms containing all three 
phases of binary coalescence, inspiral, merger and ring- 
down. 



A. Reconstruction Errors 

Errors in reconstructing these waveforms come in two 
types: errors due to SVD truncation, and errors due to 
reconstruction coefficient interpolation. The truncation 
errors have previously been shown to take the form 



Sp(x,y) 
p(x,y) 



1 N 

li=N' + l 



IM^yM 2 , (9) 



where the sum is over the basis vectors that are discarded. 
The interpolation errors have a similar form 



p(x,y) 



interp 



1 N ' 

-Y,\M»(x,y)-M' fl (x,y)\ 2 . (10) 



It should be noted that here the sum is over the basis 
vectors that are kept from the SVD. By setting the re- 
construction coefficients with \i > N f to zero, these can 
be combined into a single expression 



N 



A 4 — 1 



III. RESULTS 

We apply this procedure in two ways. In section III A 
we investigate using this approach in the context of in- 
terpolating whitened waveforms. This would be useful in 
the context of parameter estimation. Specifically, one 
could obtain reconstruction coefficients that would be 
used for constructing filter outputs associated with arbi- 
trary points in parameter space using the filter outputs 
from the SVD basis vectors. 

In section III B, we apply similar techniques to interpo- 
late raw waveforms. This is done in the context of wave- 
forms one would receive from numerical relativity simu- 
lations (i.e., time series of ^2^) = 9fh^{t) -\- id^h x {t) 



3 




FIG. 1. Reconstruction coefficients as a function of M and q 
associated with the 3 rd and 21 st basis vectors are shown in the 
left and right columns, respectively. The top row shows the 
real part of the reconstruction coefficients. The bottom row 
shows the imaginary part of the reconstruction coefficients. 



FIG. 2. The upper panel shows the fitting factor residual 
associated with using the interpolation net as templates in 
a template bank. The lower panel shows the interpolation 
mismatch for waveforms from with the parameter space. The 
interpolation mismatch is more than an order of magnitude 
smaller than the fitting factor residual. 



that are restricted to lie along lines of constant M) . This 
approach could be taken to extend numerical relativity 
waveform catalogs at greatly reduced computational cost. 



A. Whitened waveforms 

We apply this procedure to non-spinning phenomeno- 
logical inspiral-merger-ringdown (IMR) waveforms [17] 
with M E [6OM ,8OM ], q £ [1,10], whitened with 
an initial LIGO amplitude spectral density (ASD), and 
transformed to the time domain. We generate a stochas- 
tic template bank [ ] with 99% minimal match for this 
range of parameters. Since we are working with IMR 
waveforms, there is no well denned end of the waveform. 
We choose to align the whitened waveforms according to 
their peak amplitudes and compute the SVD basis vec- 
tors from these waveforms using the procedure described 
in [11]. At this intermediate stage, if we were to look at 
how the resulting reconstruction coefficients vary in pa- 
rameter space, we would see high frequency features that 
would be difficult to resolve with interpolation without a 
high density interpolation net. 

Fortunately, these features can be ameliorated by a 
complex rotation of the input waveforms, which is equiv- 
alent to a complex rotation of the reconstruction coeffi- 
cients, 

M^x,y) -»■ e-^^^M^y). (12) 

This rotation is chosen such that ^ [Mi(x, y)] = 0. Fig. 1 
shows the reconstruction coefficients associated with the 



3 rd and 21 st basis vectors after this complex rotation. 
The smoothness of these reconstruction coefficients indi- 
cates that interpolation should be possible. 

In order to perform the interpolation, waveforms from 
the (12,12) order 2D-Chebyshev net are then projected 
onto these basis vectors to obtain the interpolation co- 
efficients, as described by (7), and rotated as described 
above. 40 x 40 test waveforms from within the param- 
eter space, laid out in a grid, are used for computing 
mismatches between the interpolated waveforms, given 
by (8), and the original waveforms. Fig. 2 compares the 
fitting factor residual, which we define to be one minus 
the fitting factor, obtained from using the net waveforms 
as templates with the interpolation mismatches associ- 
ated with the test waveforms. We see that the largest 
interpolation mismatch is more than an order of magni- 
tude smaller than the fitting factor residual from the net 
waveforms. 



B. Raw waveforms 

We apply similar techniques to waveforms of a type 
that would be provided by numerical relativity simu- 
lations. Specifically, we use non-spinning phenomeno- 
logical IMR waveforms [ ] with M G [6OM ,8OM ], 
q G [1,6], multiplied by / 2 , which is equivalent to tak- 
ing two time-derivatives, and transformed to the time 
domain. We use the same alignment and rotation tech- 
niques described in section III A to prepare the wave- 
forms for interpolation. 
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FIG. 3. The upper panel shows the locations for which wave- 
forms were produced, chosen by a stochastic template place- 
ment algorithm. These waveforms were used in constructing 
the basis vectors enclosing this region of parameter space. The 
lower panel shows the interpolation mismatch for waveforms 
from within the parameter space. Waveforms interpolation 
accuracies are below a few times 10 -4 . 

To generate the basis vectors that enclose this pa- 
rameter space, we construct a stochastic template bank 
with an additional constraint. The mass ratios of the 
templates are restricted to take on values q G {qj = 
5xj + l\j G [1,6]}, where x j are the nodes associated 
with the 10th order Chebyshev polynomial. 

With the basis vectors in hand, we project the wave- 
forms from an interpolation net consisting of the (20,10) 
2D- Chebyshev nodes onto the basis vectors to obtain the 
reconstruction coefficients. These complex coefficients 
are rotated as described above, and then used to ob- 
tain the interpolation coefficients. Again, 40 x 40 test 
waveforms from within the parameter space, laid out 
in a grid, are used for computing mismatches between 
the interpolated waveforms and the original waveforms. 
We find comparable interpolation mismatches for these 
non- whitened waveforms, shown in figure 3, as for the 
whitened waveforms. 



IV. CONCLUSION 

Using the procedure described above, we have shown 
it is possible to produce gravitational waveforms for ar- 



bitrary points in parameter space by interpolating recon- 
struction coefficients from the SVD of a set of waveforms 
uniformly covering the space. 

Results have been presented for both whitened wave- 
forms, and raw waveforms. The former could be use- 
ful in the context of parameter estimation associated 
with compact binary coalescence (CBC) gravitational- 
wave (GW) signals, which frequently uses Monte Carlo 
Markov Chain methods to measure the likelihood ratio 
from many points in parameter space. This requires the 
generation of the waveforms for each point in parameter 
space and the overlap computation between the waveform 
and the data. Using the interpolated reconstruction co- 
efficients, the same computation could be approximately 
performed with generating a subset of the total wave- 
forms, reconstructing the overlap by appropriately re- 
combining the filter outputs from the basis vectors. The 
latter could be used to accurately interpolate waveforms 
that are computationally costly to produce, as is the case 
for numerical relativity waveforms. 

For future work, these techniques should be expanded 
to include additional dimensions of parameter space (e.g., 
binary object spin parameters). In addition, other inter- 
polation schemes that use equispaced or random points 
in parameter space might be found to be favorable for dif- 
ferent applications. We also note that these techniques 
could be applied to other gravitational waveforms such as 
supernova waveforms where singular value decomposition 
has also been applied [ ], or where other methods have 
been used to reduce the rank of the parameter space [26] . 
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